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Abstract 

We present a very simple and fast algorithm for the numerical solution of viscoplastic flow problems without 
prior regularisation. Compared to the widespread alternating direction method of multipliers (ADMM / 
ALC2), the new method features three key advantages: firstly, it accelerates the worst-case convergence 
rate from 0{1/Vk) to 0{l/k), where k is the iteration counter. Secondly, even for nonlinear constitutive 
models like those of Casson or Herschel-Bulkley, no nonlinear systems of equations have to be solved in the 
subproblems of the algorithm. Thirdly, there is no need to augment the Lagrangian, which eliminates the 
difficulty of choosing a penalty parameter heuristically. 

In this paper, we transform the usual velocity-based formulation of viscoplastic flow problems to a dual 
formulation in terms of the stress. For the numerical solution of this dual problem we apply FISTA, an 
accelerated first-order optimisation algorithm from the class of so-called proximal gradient methods. Finally, 
we conduct a series of numerical experiments, focussing on stationary flow in two-dimensional square cavities. 

Our results confirm that Algorithm FISTA*, the new dual-based FISTA, outperforms state-of-the-art 
algorithms such as ADMM / ALG2 by several orders of magnitude. We demonstrate how this speedup can 
be exploited to identify the free boundary between yielded and unyielded regions with previously unknown 
accuracy. Since the accelerated algorithm relies solely on Stokes-type subproblems and nonlinear function 
evaluations, existing code based on augmented Lagrangians would require only few minor adaptations to 
obtain an implementation of FISTA*. 

Keywords: fast proximal gradient methods, augmented Lagrangian methods, viscoplastic fluids, adaptive 
finite elements 
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1. Introduction 

Viscoplasticity is a wide-spread phenomenon in both natural and man-made applications. The rich 
rheology of viscoplastic fluids is encountered in geophysics, considering the examples of lava flows or lahars 
[HE]- Certain types of mineral oils, mud or slurry suspensions also exhibit viscoplastic features. In the 
consumer goods industry, toothpaste, hair gel, tomato sauce or dough serve as classical examples of such 
fluids [3]. 

The characteristic feature of a viscoplastic fluid is its ability to resist stress in the material up to a critical 
threshold, the so-called yield stress tq. This behaviour is generally due to friction-type interactions between 
the molecules or particles of the fluid. Consequently, viscoplastic fluids behave like a rigid material at small 
stress. They only start shearing like a viscous liquid if the stress exceeds the threshold posed by the yield 
stress. 
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1.1. Mathematical Models for Viscoplastic Fluid Flows 

We consider the problem of steady, creeping viscoplastic flow in a cavity, represented by the bounded 
domain fl C with (Lipschitz) boundary did = F. In practice, d G {2,3}. Our objective is to solve for 
functions u : il —?■ R"^, p : — > R and t : fl —>■ Rf^^i representing the flow velocity, pressure and deviatoric 

part of the stress, respectively. Furthermore, with the symmetric gradient operator V := (V + V^)/2, we 
denote the strain-rate tensor by 7 : SI —>■ Rf^^j which is linked to the flow velocity through the relation 
7 = Vu. 

The most common mathematical descriptions of viscoplastic behaviour are given by the Bingham [4], 
the Casson and the shear-thinning Herschel-Bulkley model [5] . With viscosity or consistency parameters 
p, K > 0 and an exponent 1 < r < 2 , they can be formulated as 




if 7 = 0 

( 1 . 1 a) 

2p7 + ro^ 
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if 7 ^ 0 . 

( 1 . 1 b) 


(Herschel-Bulkley) 




Here, | • | denotes the Frobenius norm on Rfy^. 

In what follows, we consider a non-dimensionalised formulation that has been re-scaled with respect to 
a characteristic length L and velocity U, which reduces the dimensional constitutive relations (O) to 


t\ < Bi 
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(7^ + CBi)^F 

(Casson) 

if 7 ^ 0 . 
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The Bingham number Bi := (Bingham, Casson) or Bi := (Herschel-Bulkley) quantifies the 

deviation of the viscoplastic flow from (generalised) Newtonian behaviour. 

Any of these constitutive relations, along with equations for conservation of momentum and mass, yield 
a system for the unknown flow variables. Denoting by / : H > R"^ a non-dimensionalised density of body 
forces, we have 


— Divx-|-Vp = / in H (1-3) 

divM = 0 in H. (1-4) 

To close the system, we incorporate the boundary condition 

u = md on F, (1-5) 

where Mq : Fq —> R'^ is given. We use the notation div (resp. Div) for the (rowwise) divergence operator. 
1.2. Variational Formulation 

In the following, we use boldface letters for spaces to denote d-fold Cartesian products, e.g. for a space 
A we write A := A‘^. To obtain a mathematically rigorous formulation of the viscoplastic flow problem 
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(1.2)-(1.5) in Sobolev spaces, we consider 


U 


U< 


0 * 


t/*o 

C^oo 


= W^’^(n) 

= {uG I div'it = 0} 

= {ueW^'^in) |M|r = 0} 

= {u G W^’’’(n) I divM = 0 and M|r = 0 } , 


with r = 2 in the Bingham and Casson settings. We use the dual of the latter space to fix the inhomogeneity 
/ € Uqq and we pick boundary values md G C^d, where 


17 d := <^ md G 


md • ds = 0 


Furthermore, we dehne the convex set of admissible solutions 

t^oD := { M G W^’’'(n) I divM = 0 and M|r = Md } • 

For the strain-rate and stress tensors, we will also need spaces of symmetric matrices whose entries satisfy 
an integrability condition of order r or r* , respectively, where 1/r -\- 1/r* = 1: 

By generalising the ideas of Duvaut and Lions [71 [S] and Huilgol and You [5] , we conclude that the system 


(|1.2|)-(|1.5|) is a strong formulation of the following variational inequality problem of the second kind: find 

( 1 . 6 ) 


u G UoD such that for all test velocity helds v G Uod 

a{Vu,Vv - Vu) + j{Vv) - j{Vu) > {f ,v - u)u-^,Uoo- 
This variational inequality is composed of the elliptic form a : Q x Q —>■ R, 


a( 7 , ^) := 2 / 7 : ^ dx 


0(7, ^) := 2 J 7 : ^ dx -I- 2-\/2Bi 
a I 

0(7, := 2’'“^ J |7r“^7 : ^ dec 

n 

the nonsmooth functional j : Q —R, 


7 


VWl 


: 6 da; 


(Bingham) 

(Casson) 

(Herschel-Bulkley) 


j{i) ■= Bi J | 7 |da; 
a 

and, on the right-hand side, a duality pairing between C/qo and its dual, which can be represented as 

if, V - u)u-^,Uoo = j f ■{V-U)dx 
Q 

provided that f G (fl). The colon represents the Frobenius inner product of two dx d matrices, the dot 
the scalar product of two vectors in R^^. 
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It is an important observation that for each of the three viscoplastic models, the term a{'Du,'Dv — Vu) 
possesses special structure: with the functional 6 : Q —)■ R dehned by 


Hi) ■= j l7pda; 


n 

Hi)--= + 


dx 


n 


n 


Hi) :=^ Jliidx 


(Bingham) 

(Casson) 

(Herschel-Bulkley) 


we may write 

a{'Du, Vv — Vu) = (Vu6(I?m), 'D{v — u))q, q 

as a directional derivative of 6 o I? at tt in direction v — u. Consequently, we may identify the varia¬ 


tional inequality (1.6) as a first-order optimality condition of the convex, and hence equivalent minimisation 
problem 

(VP) 


min b{Vu) + j{Vu) - (/, (Q),ir(n) + ic/oD («). 

u^U ^ ^ ^ 


where the indicator functional 


^c/oD HH 


0 


if M e t/oD 


-|-oo if M ^ t/( 


OD 


enforces the incompressibility constraint (1.4) and the Dirichlet boundary condition (1.5). 


For full details of the derivation of Problem (VP), and results regarding existence and uniqueness of 
solutions, we refer to [TOl Ch 4]. 

While the Bingham and Casson flow problems are posed in Hilbert spaces, the Herschel-Bulkley model 
demands for a mathematical treatment in more general Banach spaces. Despite extra theoretical challenges, 
a very practical consequence of this fact is that a numerical optimisation algorithm would include a series of 
nonlinear subproblems since the exponents r and r* are different and no longer equal to two. To circumvent 
these difficulties, it has become common practice in viscoplasticity (see e.g. the appendix of [9] and the 
references therein) to discretise the problem first, and then solve a discretised approximation of ( |VP[ ). This 
way, all spaces U, Q, S are replaced by finite-dimensional spaces Uh, Qh, Sh, e.g. spaces of finite elements, 
for which the Hilbert-space structure can be recovered. 

Strictly speaking, we would now have to differentiate between such finite-dimensional spaces, whenever 
we refer to the Herschel-Bulkley problem, while we concurrently work with the original function spaces 
for the Bingham and Casson problems. In an attempt to provide a clearer picture, we will not show this 
distinction in our notation and hide the subscript h for now. We shall however emphasise that in order to re¬ 
establish a mathematically rigorous formulation for Herschel-Bulkley fluids, one has to replace all variables 
with their hnite-dimensional counterparts. 


1.3. State-of-the-Art Approaches to Solving (VP) 


The community of numerical analysts in viscoplasticity can essentially be divided into two groups: some 


authors approximate the constitutive relations (1.2) with a more regular formulation, while others leave the 


genuinely nonsmooth nature of the original problem (VP) unaltered. 


For formulations of the first kind, we exemplarily mention the Bercovier-Engelman model m, the 
Papanastasiou regularisation m, bi-viscosity formulations of Tanner and Milthorpe m and De los Reyes 
and Gonzalez Andrade [Huniisiin] and the penalty approach of Glowinski and collaborators [niiiiiin]. 
It can be seen as the main advantage of these approximations that they allow for very efficient numerical 
methods of Newton-type, with a fast, locally superlinear or quadratic convergence rate. However, for many 
practical applications it can be problematic that such approximate solutions do not generally reflect the 
characteristic features of the exact solution. For instance, Moyers-Gonzalez and Frigaard [H] demonstrate 
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that even under arbitrarily small excitations, solutions to regularised models predict slow flow although 
the actual flow rate is exactly zero, see also [22l[23]. Another well-known issue arises from the difflculties 
of recovering the yielded and unyielded flow regions under a regularised problem formulation. In fact, it 
appears that convergence of these approximate plug and shear regions has not been proved yet, and such a 
statement may not even hold |23| . 


In our work, we employ the genuinely nonsmooth formulation (VP I. Classical algorithms for the numer¬ 


ical solution of this convex optimisation problem stem from the framework of augmented Lagrangians. A 
set of augmented Lagrangian methods ALG1-ALG4 was originally proposed by Fortin and Glowinski [24] . 
They are based on the idea to introduce 7 = 'Du as a constraint to the minimisation problem, the violation 
of which is then penalised with an extra quadratic term: 


min 6 ( 7 )-bj( 7 )-(/,w)i,^ 

{u,-r)eUxQ 


‘(Q),L^(n) + ''(7od('“') + f II7 “ 


subject to 7 = Du. (VP') 


The problems (VP) and (VP'„) are clearly equivalent for any p > 0. 


As another equivalent alternative to (VP' ), one frequently encounters the formulation 


, .min b(Du) + J( j) - (/, u)^, 
(u,~i')eUxQ 


‘(n),L'-(n) + ‘'Uooiu) + ^ II7 - ^“IIq 


subject to 7 = Du 


in the literature on viscoplastic flow, or the corresponding first-order optimality conditions. We refer to 
[71I1II1IIH1 for a more detailed discussion in the Bingham setting. This approach has the slight disadvantage 
that it leads to a rather complicated dual problem. As shown, e.g. in [TB], it has the form of an elliptic 
optimal control problem with pointwise inequality constraint on the control. The dual multiplier has the 
dimension of a stress and can be interpreted as a plastic contribution to the extra stress tensor. We will 
consider the formulation in (VP'^), where the Lagrange multiplier can be identified with the stress t itself. 
^^Gl, the generalised Lzawa method [IS] applied to the augmented Lagrangian that corresponds to 
attempts to find a saddle point by minimising in the variables {u, 7 ) and then taking a step along 


(vp: 


the dual gradient in order to maximise with respect to the Lagrange multiplier t. Since minimising jointly 
in u and 7 is as difficult as solving the original problem, every iteration of ALGl is very costly. 

In ALG2, the alternating direction method of multipliers (ADMM), the exact solution of each subproblem 
is waived in favour of a minimisation in u only and a subsequent minimisation in 7 only. ALG2 is related to 
the Douglas-Rachford splitting algorithm. For further details, we refer to Glowinski’s recent review [521 , the 
references therein as well as HZ). Of all augmented Lagrangian methods, ALG2 is by far the most popular 
one and can be considered as the standard genuinely nonsmooth approach to simulating viscoplastic fluid 
flows. Even though the convergence analysis of ALG2 is difficult, it is meanwhile well-known that even 
under assumptions that are too strong for Problem ( |VP[ ), only a sublinear convergence rate of 0{l/k) in 
the dual objective functional can be guaranteed [551 Thm 1]. This corresponds to a convergence rate of only 
0 (l/-\/fc) for the primal iterates, which represent the velocity and strain rate. 

ALG3 is a counterpart of the Peaceman-Rachford splitting method. Also known as alternating minimi¬ 
sation algorithm (AMA), it expands on ALG2 by adding an update of the dual variable r not just after the 
minimisation in 7 , but also after the minimisation in u. ALG4 was designed as an analogue of the 0-method 
[531152], but appears to be rarely applied. 


Gonvex optimisation problems of a very similar structure to (VP) are also encountered in a variety of 


other disciplines, such as signal and image processing [301E], machine learning, statistics or mathematical 
finance. The field of nonsmooth convex optimisation consequently receives considerable attention and many 
alternative numerical methods have been derived in past years. To us, the class of so-called proximal gradient 


algorithms appears particularly well-suited for the solution of Problem (VP), as these methods can readily 


exploit the structure of composite convex objectives: terms which are smooth in the sense that they possess 
a Lipschitz-continuous gradient are linearised, while nonsmooth terms are left unchanged. 

In recent years, interest in proximal algorithms has reached an unprecedented extent with the re-discovery 
and further development of fast or accelerated methods. Building upon ideas of Nesterov [35[|33], Beck 
and Teboulle [M] constructed a fast iterative shrinkage-thresholding algorithm (FISTA), an efficient fast 
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variant of the basic proximal gradient method. At every iteration, FISTA extrapolates from the current and 
immediate past iterate to accelerate the convergence of the method. While the conventional, unaccelerated 
method decreases the value of the objective at a rate of order 0 (l/fc), where k is the iteration counter, 
FISTA achieves 0{l/k'^) at negligible extra computational cost. In particular, it requires no accumulation 
of the iteration history, apart from the last iterate. 

Very recently. Beck and Teboulle considered abstract composite problems of the form 


min/(M) + g{'Du) 

U 

with a strongly convex functional / and a general convex functional g. In [35] (see also |36j 1. they apply 
FISTA to the dual problem to achieve 0(I/fc)-convergence of the primal sequence, in contrast to 0(I/-\/fc) 
without acceleration. Given that problem (IVPl) is of the similar form 


min fiT>u) + g{u) 


again with a strongly convex / and convex g, our approach exhibits some analogies to |35j . The similarity 
between our and their problem allows us to apply some of their ideas in our work, in particular the concept 
of applying FISTA to the dual problem. However, there are also two important differences: firstly, the 
different structure of the problem gives rise to some additional difficulties and requires further technical 
assumptions in order to guarantee the applicability and convergence of FISTA. Secondly, Beck and Teboulle 
restrict themselves to problems in R", while we consider a general function space setting. 

Even though our approach is applicable to very general composite convex problems, we will only focus 


on the application of the methodology to Problem (VP). The more abstract mathematical details can be 
found in [ 10 ]. 

We shall point out that acceleration techniques have also been studied for the augmented Lagrangian 
methods ALG2 and ALG3, see [551 157] . However, quantitative rates of convergence are very difficult to 
prove for these methods. At this stage, it appears that there are only isolated results available that rely on 
fairly restrictive assumptions or unconventional measures of convergence (28l|38l|39l|40l|4T]. 

For an extensive review of state-of-the-art gradient-based optimisation methods for abstract composite 
convex problems, we refer to [42] . 

Outline. In Section]^ we introduce a novel dual formulation of Problem ( |VP[ ). For its numerical solution, we 
present the accelerated dual proximal gradient method FISTA* in Section [^ along with some key properties. 
Finally, we conclude with numerical results and a discussion in Section]^ 


2. Dual Formulation 

2.1. Derivation of the Dual Problem 

We consider Problem ( VPq ), the split formulation of the viscoplastic flow problem with no penalty term 
added: 


min 6 ( 7 )-h j(7) - +(-( 7 odN subject to 7 = (VPq) 

(u,-V)GC/xQ ^ 

By introducing a Lagrange multiplier (dual variable), which has a physical interpretation as an admissible 
stress tensor r [3], we can re-write this problem equivalently as 

^ -(/,«)z,'-*(n),i7n)+7-2 ^«)q*,q- 

tGS iu,-y)eUxQ \ \ I ^ 

We now solve the inner minimisation problem to eliminate the primal variables u and 7 . Re-arranging yields 


max 

tGS 


min 

\-yeQ 


{ b{i)+jii) - ( t , 7 ) q*,q } 


mm 

uGU 


I -(/) “)z,'-*(n),i’'(n) + ''C/od(“) + • 


( 2 . 1 ) 


6 




For the second minimisation problem in (2.1) we infer 


mini -(/,«) 

u£U K. 




+ LUouiu) + {t,Vu)q^ q I 


0 


if X G C 


—oo otherwise. 




with the set 


C := I T G S' {t,Vu)q^ q - (/,M)i,r-*(n)_i.(n) = 0 for all m G C/qo } ■ 


We point out that this condition on r is a very compact variational formulation of (1.31. This becomes 
more obvious if we introduce the pressure p as a Lagrange multiplier for the constraint div it = 0, which is 
implicitly contained in the definition of the space Uoq. Then the set C contains all t G S, for which there 
exist p & U’ (n) such that 

{t,Vu)q^ q - (p,divit)i..(n),L.-(n) = (/,-“) 
for all test functions it G t/*o- 

Let us now turn to the first minimisation problem in (|2.1|). Separate calculations for the three different 


constitutive models reveal the following explicit result [TOl Ch 6]: 

mm{ &(7) + j'(7) “ ('r,7)Q*,Q } = --S(t) 


with 


F{t) = 


|r| — Bi)l da; 


(Bingham) 
a 

J + dx (Casson) 


( 2 . 2 ) 


^ / (kl-Bi); dx 


(Herschel-Bulkley). 


Here, we employed the notation (•)+ = max{ 0, • }. 


In summary, (2.1) leads to the following dual formulation of the viscoplastic flow problem: 


minF'(x) + lc{t) 
-res 


(VP* 


or simply 


minF'(T). 

-rec 


is a composite convex optimisation problem, where both terms 


2.2. Properties of the Dual Problem 

The dual viscoplastic flow problem (VP* 
possess special properties. 

For all three constitutive models, the functional F is (Frechet-)difFerentiable. After short calculations 
we find that its gradient with respect to r is given by 


VF(t) 


^^(M-Bi) ; 


1 
2 
1 

12 


+ I- 

\ 2 
+ , 

r* — 1 

|t| 


(Gh-v®))^ 

(ii-i - Bi); 


(Bingham) 

(Casson) 

(Herschel-Bulkley). 


(2.3) 
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In fact, it holds that WF{t) = 7 : by re-arranging this equation for the stress r, one recovers the constitutive 


models of (1.2). 


Additionally, for Bingham and Casson fluids, VF is globally Lipschitz-continuous, where the smallest 
possible Lipschitz constant is L = 1/2. This could either be verified directly, but it follows more easily from 
classical relations between primal and dual problems in convex optimisation nni Lemma 3.6]. 

Note that since 1 < r < 2 for shear-thinning Herschel-Bulkley fluids, the dual exponent r* > 2. Therefore, 
in this case VF is only Lipschitz-continuous on bounded subsets of the stress space S. However, as long as 
an algorithm for the solution of Problem (VP*) does not diverge, all stress iterates will remain bounded and 


we can still find a Lipschitz constant L > 0 for the set of all iterates. In contrast to the Bingham and Casson 
settings, it appears impossible in practice to determine this constant a priori. Instead, for Herschel-Bulkley 
flow problems we will have to compute estimates for L numerically. 

The second term in (VP*) is an indicator functional, which imposes the momentum equation in form of 
the constraint t € C. In convex optimisation, indicator functionals of convex sets serve as a prototypical 
example for functionals of simple structure. A functional G is said to be simple if its so-called Moreau 
proximal map |43j . 

proxG(T) := argminj G(cr) + \\W - t||| | , (2.4) 

can be evaluated in a computationally efficient manner. If G is an indicator function of the set G, then its 
proximal map equates to a projection of r onto G. 


We are now in the position to proceed with the numerical solution of (VP 


3. Accelerated Dual Proximal Gradient Method 


Our objective for this section is to present algorithmic approaches to solving the primal viscoplastic flow 
problem (VP I through its dual (VP*). First, we briefly review the fundamentals of FISTA and then apply 
this method to (VP*), in order to derive the basic algorithm FISTA* for the solution of (VP). After a 


comparison with the alternating direction method of multipliers (ADMM / ALG2), we present our results 
on convergence rates. Finally, we comment on some aspects of the implementation. 


3.1. Proximal Gradient Methods 

The fast iterative shrinkage-thresholding algorithm (FISTA) belongs to the class of proximal gradient 
methods. Its name stems from the fact that for certain convex functionals G, the proximal map ( 2.4[ ) 
becomes a shrinkage-thresholding operator. However, it has meanwhile become common practice [84] to 
slightly misuse the term iterative shrinkage-thresholding algorithm even if the proximal map of G is of 
different form, as is the case when G = be actually yields an iterative projection algorithm. 

The most basic proximal gradient method, or ISTA, for the solution of the composite convex problem 

minF(T) -|- G(t), 

-tGS 

where F has a Lipschitz-continuous gradient with Lipschitz constant F, reads as follows |34j : 

Algorithm 3.1 (ISTA / Proximal Gradient Method). Input: S S 

Initialisation: k = 1 

(ISTA.l) Set = L or find a valid > 0 by backtracking and evaluate 


(ISTA.2) 
(ISTA.3) 


If the algorithm has converged Then Return and Stop. 
Set fc •(— fc -I- 1 and Go To (ISTA.l)"] 
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A valid estimate of the Lipschitz constant L means that and the corresponding iterate 
must satisfy the descent criterion 




L(fe) 


^(fc) _ ^(fe-i) 


(3.1) 


which holds in particular for all > L [Ml Lemma 2.1]. To find a valid parameter in an algorithmic 
fashion, one would start from a trial value L and increase this value by a certain percentage until ([33 is 
satisfied. We will present the full details below in Algorithm |3.4| 

The following convergence result for ISTA is established in |M], Remark 2.1 and Theorem 3.1: 

Theorem 3.2 (Convergence Rate of ISTA). Let the sequence (t^^))^ be generated by ISTA, where the 
sequence is non-decreasing and bounded. If t G argmin.,.gg { F{t) + G{t) } is any solution of the 

minimisation problem, then 


{F{t) + G{t)) - (Fir) + G{t)) < 0{l/k). 

He and Yuan |39j prove a corresponding 0(l//c)-convergence result for the alternating direction method 
of multipliers (ADMM / ALG2), but it appears that their stronger assumptions on the problem compared 
to Theorem |3.2| cannot generally be relaxed. 

It turns out that this convergence rate of ISTA and ADMM / ALG2 for the class of problems with smooth 
convex F and convex G is not optimal in the sense of complexity theory (cf [HI pp 4-7]). In fact, there 
exist algorithms which converge like 0(l/fc^), while it can be shown |Mj that no higher rate is achievable 
for this entire class of optimisation problems. The fast iterative shrinkage-thresholding algorithm FISTA 
constitutes an example of a method which is optimal in this sense. Such accelerated or inertial algorithms 
have only recently attracted greater attention in convex optimisation. 

Compared to ISTA, FISTA comprises an additional extrapolation step. It combines the last two iterates 
and in just the right amount to compute a so-called leading point For the next iteration, 

the functionals F and G and the gradient VA are then evaluated at this leading point instead of 
to find a shortcut to the minimum. 

Algorithm 3.3 (FISTA / Accelerated Proximal Gradient Method). Input: G S 

Initialisation: k = 1, t^^'^ = 1, 

(FISTA.1) Set = L or find a valid > 0 by backtracking and evaluate 

xW = prox^G - ;^VF(tW)) 


(FISTA.2) If the algorithm has converged Then Return and Stop. 
(FISTA.3) Gompute 

^ 1 + Vl + 4t(fc)2 
2 

to update the leading point 




t(k) _ 1 



(FISTA.4) Setk^k + l and Go To (FISTA.l) 


Since FISTA introduces the additional variable x^^\ its memory footprint is larger than the one of ISTA. 
However, the computational cost of the extrapolation step (FISTA. 3)j which does not appear in ISTA, is 
only marginal: it solely requires the evaluation of a linear combination. 
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Similar to (3.1), if no Lipschitz constant L is known, then the estimate must at least satisfy the 

lW 


descent criterion 


' ' s* ,s ^ 


T-ik) _ ~{k) 


(3.2) 


We recall that the Bingham and Casson models allow us to set L^^'> = L = 1/2. For Herschel-Bulkley fluids, 
the following backtracking strategy [31] can be employed in step (FISTA.l) 


Algorithm 3.4 (Backtracking). Input: a trial value for > 0, a magnifying factor rj > 1 

(BT.l) Evaluate 


= prox 1 n 

l(>=) 


T-(fc) _ 




VF(t 


(k)- 


(BT.2) If (3.2) holds Then Return and Stop. 


(BT.3) Set ^ and Go To (BT.l) 


For both fixed and variable Beck and Teboulle derive the following result [Mj (Remark 2.1 and 

Theorem 4.4): 

Theorem 3.5 (Convergence Rate of FISTA). Let the sequence ('r^^))^ he generated by FISTA, where the 
sequence is non-decreasing and bounded. If t G argmin.,.gg { F{t) + G{t) } is any solution of the 

minimisation problem, then 

{F{t) + G{t)) - {F{t) + G(t)) < Oil/e). 

3.2. The Accelerated Dual Proximal Gradient Method FISTA* 

Let us turn towards the application of FISTA to the solution of viscoplastic flow problems. We will now 
derive how the proximal map in step (FISTA.l) can be evaluated in practice and how the primal variables 
u and 7 , representing the flow velocity and strain rate, can be recovered. 

Lemma 3.6. The assignment = prox^Q(T*'^^ — WF{t^^'^)/L) is equivalent to the operations 

7 ^^^ = arg min | 6 ( 7 ) + j( 7 ) - (r 7 ) j 
-ygQ I ' ' Q*.Q J 


;(fc) _ 


u£U 


= arg min <{ -(/, M)i,r-* + iuoui^) + ^ 


( 


Vu- ( 7 ^''^ -Lf 


“")ii:} 


r(C = + 1 {vu^k) _ . 

Proof. We refer to [H Ch 3]. 


(3.3) 

(3.4) 

(3.5) 

□ 


In (3.3) and (3.4), the quantities 7 ^ ^ and represent a strain rate and a velocity, respectively. They 
are, however, evaluated based on the extrapolated leading point not the actual stress iterate The 
convergence result for FISTA is valid for this sequence (x*^^))fc only, while such a convergence result may 
not hold for the auxiliary sequence 

We suggest three definitions for a primal sequence [u^k) of velocity and strain-rate fields: 

• The first idea is to solve the problems in (3.3) and (|3.4[) once again, this time with instead of 


A^^k) ^ arg min I 6 ( 7 ) -f j( 7 ) - (T^k)^A^\ \ 

-VgQ I ' ' Q*,Q J 

-(f,u}Lr-*^n),L’-(n)+kUon(u) + ^ Du - (y^k) _ 


u^k) = arg min 
uGU 


(3.6a) 
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• The second idea is to obtain from 7 ^^^ by solving the equation = 7 ^^^ in a least-squares 

sense: 


= argminj 6(7) -f j(7) - 1 

-VgQ L ' ' Q*,Q j 


= argmin 
uGI/qd 


Vu — 7 


(fe) 


(3.6b) 


• The third idea is to simply set 


\^{k) 


(3.6c) 


Combining Algorithm |3.3[ Lemma |3.6| and these definitions of the primal sequence, we obtain the accel¬ 
erated dual gradient method FISTA* for solving (VP). 

Algorithm 3.7 (FISTA* / Accelerated Dual Proximal Gradient Method). Input: G S, gradTol > 0 

Initialisation: fc = 1, = 1, 

(FISTA*.1) Set = L or find a valid > 0 by backtracking and evaluate 


c(fe) 


-yeQ 


7 ''"'= argmin <; b(j) + j(j) - 


Q*,Q 


(FISTA*.2) Solve 


= argmin 
u£U 


I -if, {n),L^{n) + kuon («) + - (t ||^ | • 


(FISTA*.3) Update 




(FISTA*. 4) Obtain and from (3.6). 


(FISTA*.5) If —-y(k) £ gradTol Then Return and and Stop. 

Q 


(FISTA*.6) Compute 

to update the leading point 


^(^+1) _ 


1 -F vTT4t(U2 


dfc+i) = T-(fe) , 1 /, 

tik+i) V 


.^(fc) _ .^(fc-i) 


)■ 


(FISTA*.7) Setk^k + 1 and Go To (FISTA*. 1) 


Let us briefly analyse the minimisation problems in (FISTA*.1) and (FISTA*.2) from properties of 
convex conjugates it follows immediately [TOl Ch 2] that [(FISTA*. 1) explicitly reads 

= VF(r(''^) 


II 


























with VF from (|2.3|). 


and € F’’ (O) such that 


(FISTA*.2) is a variational formulation of the following Stokes problem: find S C/qd 


1 

'Zw 


DivF'ii*'^^ + = f — Div ^ ~ 


div = 0 . 


Analogously, the least-squares problem in (3.6b) leads to the Stokes problem 

- Div = -DivV'^) 

div = 0 . 


3.3. Comparison with Classical Algorithms 

Since conventional algorithms do not include the extrapolation step to determine a leading point, let us 
consider the non-inertial counterpart of the dual FISTA method, i.e. the dual ISTA method. In the absence 
of a leading point, defining a separate primal sequence becomes redundant. 

Algorithm 3.8 (ISTA* / Dual Proximal Gradient Method). Input: € S, gradTol > 0 

Initialisation: k = 1 

(ISTA*.l) Set = L or find a valid > 0 by backtracking and evaluate 

= argminj b{j)+j{j) - 

■yeQ I ' ' Q*,Q 


(ISTA*.2) Solve 




argmin{-(/,M)i.*(n),i-(n)+i; 7 oD(«) + ^|- 


(ISTA*.3) Update 




Li^) 




(ISTA*.4) If < gradTol Then Return and and Stop. 


(ISTA*. 5) Setk^k+1 and Co To (ISTA*.l) 


The alternating direction method of multipliers applied to the formulation (VP), 


reads as follows: 


Algorithm 3.9 (ALG2 / Alternating Direction Method of Multipliers). Input: 7 ^°^ € Q, S S', a 
positive sequence of step sizes, gradTol > 0 

Initialisation: k = 1 


(ALG2.1) Solve 


ti(G =argmin| -(/,M)i,.--(o)^i 7 n) +'-c/od(“) + | 


Q 


(ALG2.2) Solve 


7 ^^) = argmin I 6 ( 7 ) + j( 7 ) — +- 1 . 

yeQ I ' ' Q*,Q 2 q J 
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(ALG2.3) Update 








7 


(fc) 


Q 


< gradTol Then Return and and Stop. 


(ALG2.4) If - 7^'=) 

(ALG2.5) Set fc ^ fc + 1 and Go 7d |(ALG2.l)1 

We observe a number of similarities and differences: 

• The steps (ALG2.1) and (ISTA*.2) are identical, provided that g = 


The steps (ALG2.2) and (ISTA*.l) are identical, provided that g = 0. It can be seen as a ma¬ 


jor disadvantage of the alternating direction method of multipliers with g > 0 that for Gasson and 


Herschel-Bulkley fluids, the solution of (ALG2.2) cannot be obtained explicitly. Instead, an iterative 


scheme such as Newton’s method is required in every iteration of ALG2 to compute 7 ^^\ In (ISTA*.l) 


ISTA* only evaluates a nonlinear function. The same holds true for (FISTA*.l) 


The updates (ALG2.3j|(ALG2.'5) and (ISTA*.3)|(ISTA*T5) are identical, provided that 


• The parameters g and in ALG2 have to be chosen heuristically, whereas a globally optimal value 
for can be calculated analytically, or estimated by backtracking. 

In conclusion, ISTA* solves either the same or simpler subproblems than ALG2 in every iteration. 
FISTA* additionally incorporates an extrapolation step to achieve a higher rate of convergence. We will 
quantify these rates for each algorithm in the next subsection. 


3 . 4 . Convergence of FISTA* 


In step (FISTA*.4) (3.6) provides three alternatives for defining the iterates and FISTA* with 


(3.6a) or (3.6b) approximately doubles the computational cost per iteration compared to ISTA*: the lion 


share of computing power is required for the execution of (FISTA*.2) and, to a far lesser extent, (FISTA*. 1) 


The additional step (FISTA*.4) demands for a solution of these two problems with different data a second 


time. In contrast, if the primal sequence (3.6c) based on the leading point is chosen, then (FISTA*.4) 
becomes trivial and no extra problems have to be solved. 


Unfortunately, convergence of the sequence (3.6c) cannot be guaranteed. At this stage, it is an open 


problem whether the convergence of this sequence might actually fail, whether it generally converges even 
though no proof appears to be available, or if additional assumptions on the problem are required to ensure 
convergence (cf [101 Ch 3]). 

For the other two alternatives, we have the following result: 


Theorem 3.10 (Gonvergence of FISTA* and ISTA*). Let u he the exact solution of Problem ( |VP| ). 

If the sequences ( 7 ^^^)fc and are generated by FISTA* with either (|3.6a[) or (|3.6b[), then 


7 *-^) — Vu 


- u 


Q 


= 0(I/fc) 


^=0(I/fc). 

If the sequences (7^*^^)fc and are generated 


by ISTA*, then 


-yi^) _ 


— U 


= 0{1/Vk) 
= 0{1/Vk). 


Proof. We refer to [THl Gh 3]. 


□ 
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Let us mention once again that under additional regularity assumptions on the problem, the alternating 
direction method of multipliers ALG2 with a suitable choice of parameters exhibits the same convergence 
rate as ISTA. 

The sequence of dual multipliers is bounded and therefore possesses a weakly convergent sub¬ 
sequence. If weak convergence of the entire sequence is desired, then this can be achieved by a 

modification of the extrapolation sequence in FISTA* [i5] . 


3.5. Computational Techniques 

Solution of the Stokes Problems. For solving the Stokes problems that occur in each method, we apply the 
classical preconditioned conjugate gradient Uzawa (PCGU) method of Cahouet and Ghabard [46l pp 892- 
893]. Glowinski [THl Sec 20-22] motivates this method by successively improving on the very basic Uzawa 
(i.e. dual gradient) method for the Stokes problem. These step-by-step improvements on the speed of 
convergence of the algorithm are achieved by 

• using conjugate gradients instead of steepest descent. 

• performing an exact line search at every iteration instead of using the global Lipschitz constant of 
the dual gradient as a worst-case estimate. Since the problem is quadratic, the exact step size is 
straightforward to calculate. 

• preconditioning the problem. While this is essential for the instationary generalisation of the Stokes 
problem, it is not beneficial for the stationary case that we consider here. With no time dependence, 
the suggested preconditioner would simply degenerate to a multiple of the identity. Looking at the 
discrete problem, one could interpret the occurrence of a mass matrix as preconditioning, though. 


We terminate the PGGU algorithm as soon as the dual gradient is sufficiently small, as measured by a 
positive constant stokesTol: 


divM('=) 


< stokesTol. 

L’-(n) 


Our convergence analysis of the dual proximal gradient methods assumes that the proximal map is 
evaluated exactly. Therefore, we should at least set stokesTol ^ gradTol. This poses no major obstacle, 
since the PCGU algorithm achieves a linear rate of convergence for the Stokes problem, compared to the 
sublinear rate of the outer optimisation loop. 


Adaptive Re-Starting. In general, FISTA is not a monotone method in the sense that the value of the dual 
objective F-\-G may also increase from one iteration to another. Analogously, the primal sequence generated 
by Algorithm FISTA* may temporarily digress from the solution {u,j). This is a well-known property of 
accelerated gradient schemes and can be interpreted as excessive momentum from past iterations, that causes 
the sequence to overshoot the minimiser and to converge in a spiralling motion. 

In [3^, Beck and Teboulle propose a simple modification of FISTA, called MFISTA, which guarantees 
monotonicity. In contrast to the basic FISTA scheme, it requires the functional F -|- G to be evaluated at 
every iteration. 

O’Donoghue and Candes |47] suggest to adaptively re-start the algorithm once an increase in the objective 
is detected in order to preserve monotonicity and to discard accumulated momentum of the iteration. Rather 
than observing the functional values, they showed a re-starting criterion based solely on the dual gradient 
to be similarly effective. According to this gradient scheme, Algorithm FISTA* is re-started whenever 

^(fc) _r('=-i)\ <0. (3.7) 

\ / s*,s 


In that case, would be an ascent direction for the dual functional at The authors of m 

point out that this scheme unites the benefits of increased numerical stability near the optimum on the one 
hand and, on the other hand, no extra computational expenditure: all quantities in (3.7) have already been 
computed previously. 
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By allowing for re-starts, the worst-case convergence rate decreases from 0{l/k) to 0{l/y/k). This is 
a consequence of the fact that the first step after (re-)starting is equivalent to a step in the unaccelerated 
dual proximal gradient method ISTA* and we have 


Therefore, 



T-(fc-l) 




1 




2 

Q 


^0. 


and this first step is generally accepted. 

Surely, the rationale behind re-starting schemes is that re-initialisations only occur as isolated events, 
not after every single iteration. This way, the convergence rate would remain close to 0{l/k). A re-start 
would ideally result in a shortcut towards the solution and thus decrease the error more efficiently than 
continued iterations with full momentum. For a numerical study of this effect, we refer to m- 


Discretisation. We assume that is a polygonal domain we let 7^ be a regular triangulation on n. We 
construct a finer triangulation 7^/2 by connecting the edge midpoints in each triangle T € Th- 

Following Glowinski [HI p 303], we apply the Pi-iso-P 2 /Pi element of Bercovier and Pironneau [H] for 
discretising the velocity and pressure, respectively. Pfc denotes the space of polynomials in two variables of 
degree at most k. We approximate the strain rate and stress with piecewise constant elements on the finer 
mesh. 

Overall, we consider the following finite-element spaces: 


Ph 

Qh 

Sh 


{ M?, e 0 (0)^ I M/ijr = 0 and € P?, VT e 77,/2 } 
{Ph(^c{n)\ph\T(^Pi, yPeTh} 

{7, e 1 e p3, vrer^/2} 




tDt&K yp&Th 


hl2 


}■ 


With MD,/t G C(r) n Ud serving as an approximation to uu that is linear on the triangle edges of 77/2jwe 
also introduce the convex set 


UtUM ■= {uh & 0(0)^ I = uu,h and Uh.\T € P^, VT G 77/2 } ■ 

It is well-known that these spaces are LBB-stable [IS], i.e. the Stokes problems discretised with these 
finite elements possess a unique solution. 

4. Numerical Results 

In this section, we will conduct some numerical experiments. Our objective is twofold: firstly, we inves¬ 
tigate the effect of the acceleration to demonstrate the improved rate of convergence of the accelerated dual 
proximal gradient method FISTA* compared to ISTA* and in particular the alternating direction method 
of multipliers ALG2, the current benchmark for solving viscoplastic flow problems with no regularisation. 
Secondly, we wish to verify that FISTA* computes accurate approximations, in particular with respect to 
predicting yielded and unyielded regions in the flow. A simple scheme for adaptive mesh refinements will 
help us to obtain high-fidelity approximations with great efficiency. 

We leave a few parameters unaltered for our simulations: for the dual gradient methods, we set = 
1/2. This Lipschitz parameter is kept constant for Bingham and Gasson flow simulations, while we choose 
a magnifying factor of ry = 1.1 for the backtracking procedure in the Herschel-Bulkley setting. For ALG2, 
we use the corresponding value g = = 2 for all k. We initialise the PCGU algorithm for solving each 

Stokes problem with the converged solution for the pressure of the previous run. All other initial guesses 
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Figure 1: Convergence history for Bingham flow in a force-driven square reservoir. 


for the required variables shall be 0 in each example. In the stopping criteria we set gradTol = 10 ® and 
stokesTol = 10“^^. We run our programs in MATLAB R2013a 64-bit on a laptop with Intel® Core^^i? CPU 
4x2.50 GHz. 

For our numerical experiments, we consider two different flow problems in a square reservoir: a force- 
driven and a wall-driven flow. In both cases, we dehne ft := ]0,1[^. To mesh the geometry, we proceed as 
follows: first, we generate a uniform grid of l//i x 1/h squares. Then we divide each square diagonally into 
four congruent triangles, the collection of which defines the coarse pressure grid 7^. 

In all our experiments with FISTA*, the primal sequence (3.6c I based on the leading point turns out to 
converge. Hence, there is no need to fall back to any of the computationally expensive alternatives (3.6a I 
or (l3^. 


4.I. Force-Driven Cavity 

In [T7] , De los Reyes and Gonzalez Andrade simulate the flow of a Bingham fluid, which is driven by the 
force 

f{xi,X2) := 300(a:2 - 0.5,0.5 - 

with Bingham number Hi = lOy/2 and homogeneous boundary conditions. We use the same parameters, 
and additionally consider the corresponding Casson and Herschel-Bulkley flow problems (r = 1.5). Our 
simulations are carried out on the grid with h = 1/32. 


Convergence Rates. We carry out 5,000 iterations of FISTA* to approximate the exact solution of (VP I 


UK, u 


^5,000)^ This allows us to report very accurate estimates of the true error 


\u 


(k) 


-u\\u. 


In Figure we compare the convergence of ADMM / ALG2, ISTA* and FISTA* for the Bingham 
flow problem. Additionally, we show how the convergence of FISTA* is affected when the criterion (3.71 is 
monitored to trigger adaptive re-starts of the method. 

In this setting, the accelerated methods achieve a convergence rate of order 0(l/k) within only few 
iterations. ADMM and ISTA* also exhibit their worst-case convergence rate of 0(l/-\/fc) after a start-up 
phase. FISTA* with adaptive re-starting discards its previous momentum after the iterations k = 144 and 
k = 351. Although the descent is qualitatively more consistent, the method is not strictly monotonous and 
in absolute numbers, the errors lie above those that can be achieved without re-starting. 

It is a very important observation that replacing iteration numbers with computing time on the horizontal 
axis does not visually affect the graphs of the error curves. This confirms the theoretical expectations that 
the cost of one iteration in ADMM, ISTA* or FISTA* is virtually identical. However, after 1,000 iterations, 
FISTA* has computed an approximation that is about two orders of magnitudes more accurate than the 
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error error 




^ computing time (s) 

Figure 2: Convergence history for Casson flow in a force-driven square reservoir. 




Figure 3: Convergence history for Herschel-Bulkley flow in a force-driven square reservoir (r = 1.5). 
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estimates returned by ADMM and ISTA*. Hence, even though FISTA* converges in a non-monotone fashion, 
it is hard to think of any practical disadvantages of this effect. 

Moving on to the results for the Casson fluid in Figure]^ fast optimisation algorithms show convergence 
of order 0(l//c^). Surprisingly, even ISTA* appears to attain the same rate asymptotically. In contrast, the 
results suggest that re-starts of FISTA* after the iterations with indices k = 313 and k = 919 decelerate the 
convergence rate down to 0{l/k). Nevertheless, the re-starting scheme effectively improves the monotonicity 
of the descent. After 1,000 iterations, FISTA* has arrived at a solution about ten times more accurate than 
the 1,000'^*' iterate of ISTA*. ADMM is uneconomical for solving the Casson flow problem as it requires an 
accurate solution of a nonlinear problem at every iteration. 

In Figure we present the results for the simulations of a shear-thinning Herschel-Bulkley fluid. No 
re-starts of FISTA* occurred during the first 1,000 iterations. The results are otherwise similar to the Casson 
flow problem. 

Yielded and Unyielded Flow Regions. Let us now compare the geometry of the stagnant zones as they are 
predicted by the genuinely nonsmooth methods of this article and the regularised approach of De los Reyes 
and Gonzlez Andrade in m- Even though the authors study a time-dependent Bingham flow problem, 
they report that a steady state is quickly attained. This allows us to compare their results under the 
quasi-stationary regime with ours. 

To visualise the regions of yielded and unyielded flow, it is common to plot areas where \Th\ < Bi and 
\Th\ > Bi in two different colours. However, in its basic form, this approach does not normally provide 
satisfactory results, as numerical errors near the interface become visible in the form of colourful noise. A 
traditional remedy [El |49l [50] exploits the fact that the solution tends to be more regular and converges 
a lot faster in yielded flow regions and therefore the shape of sets with \Th\ < (1 + £)Bi for a positive 
correction factor e is typically far smoother. Nevertheless, the main drawback of this postprocessing step 
is that it introduces a systematic error by overestimating the actual unyielded regions. The extra effort we 
have invested in the numerical solution by means of a genuinely nonsmooth method would partially be futile 
if we still apply some form of smoothing to the results. 

We therefore propose another means of visualising the yield surfaces: in the right-hand diagrams of 
Figures[4][^ we plot the magnitude of the extra stress tensor |x;j| in a window of ±0.1% around the Bingham 
number Bi. Stress magnitudes below the critical value appear in grey-green, yielded regions are displayed as 
blue-white. Therefore, the interface between yielded and unyielded regions, as predicted by the numerical 
solution, lies at the sharp transition from blue to green. Meanwhile, since the classification into yielded 
and unyielded is least reliable near \Th\ = Bi due to numerical errors, the width of blue and green shaded 
areas serves as an indicator of uncertainty in the identification of flow regions. As for the correction factor 
£, there is of course some arbitrariness in choosing the width of the interval around the yield stress, which 
defines the span of the colour gradients. However, this width never introduces any systematic errors into 
the visualisation, as the discontinuity in our colourbar always occurs exactly at the value Bi. 

Upon comparing the results in Figure|^with |17j . we observe major differences. The central solid region 
of approximately square shape deviates significantly from the cross-like or ±-shape computed by De los 
Reyes and Gonzalez Andrade. In our visualisation, it appears that the stress magnitude lies well below the 
yield stress, which indicates that those results should be reliable. 

In the corresponding graphs for Casson flow (Figure]^ and Herschel-Bulkely flow (Figure]^ we recognise 
a similar pattern of a rounded cross-like structure in the stress. Still, our numerical results suggest that the 
actual unyielded region extends beyond the grey cross and in fact the flow stagnates in the entire region 
that is shaped like a rounded square. 

On another hand, the flow field in the Bingham case computed by FISTA* agrees both qualitatively and 
quantitatively with the semismooth Newton method in m- This observation confirms that regularisation 
techniques are most appropriate if a simulation serves the sole purpose of finding an accurate approximation 
to the velocity field, but not necessarily of reflecting the exact sparsity pattern of the strain rate. 
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Figure 4: Flow velocity and plug zones for rotational Bingham flow in a square reservoir {h = 1/32). 



Figure 5: Flow velocity and plug zones for Casson flow. 
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Bi= 2 


Bi= 20 


Bi= 200 



iterations 


iterations 


iterations 


Figure 7: Iterations until convergence for different mesh sizes (of the coarse pressure grid) and yield stress parameters. ALG2 
(blue) and ISTA* (orange) failed to converge within 5,000 iterations for Bi = 200, therefore only the results for FISTA* (yellow) 
and FISTA* with re-starting (purple) are shown here. 


4-2. Lid-Driven Cavity 

We now move on to the popular benchmark problem of viscoplastic flow inside a lid-driven square cavity. 
The fluid motion is only driven by a moving wall, i.e. we have f = 0 and we define 

/ (1,0)^ if X2 = 1 
1(0,0)^ otherwise 

We point out that due to the discontinuities in the top corners, this choice violates the assumption md G tdo- 
The lid-driven cavity problem is classically studied to assess the performance of numerical methods under 
the presence of singularities. 


Iterations and Computing Times. For a range of different grid sizes and values of the Bingham number, 
we now compare how many iterations and how much time each of the four algorithms ADMM / ALG2, 
ISTA*, FISTA* and FISTA* with re-starting requires to compute a solution of the prescribed accuracy 
gradTol = 10““^. For these simulations, we focus on the problem of Bingham flow. 

It turns out again that by incorporating an acceleration scheme, the number of iteration is reduced 
significantly in every single case. For the largest value of the Bingham number considered here, Bi = 200, 
ALG2 and ISTA* were still far from an optimal solution even after 5,000 iterations, which we consider as 
failed to converge. It is well possible that for different algorithm parameters ALG2 would have terminated 
within these 5,000 iterations. However, such parameters would have to be found by trial and error, which 
is clearly unpractical. For the sake of comparability with ISTA* and FISTA*, where such heuristics are not 
necessary, we limit our presentation to the setting g = s^^'> = = 2. 

For the other cases, where we have data for all algorithms available, the dual FISTA method requires 
83% fewer iterations and 79% less computing time than the alternating direction method of multipliers. The 
reduction in iteration numbers and GPU times for the re-starting adaptation are 83% and 78%, respectively. 
As can be seen from Figures and re-starting is worthwhile in certain examples, while it is the opposite 
in others. 
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Figure 8: Computing times corresponding to the test runs in Figure 


J^.S. Yielded and Unyielded Flow Regions 

As noted by Yu and Wachs [15] , the precise resolution of yielded and unyielded regions becomes partic¬ 
ularly challenging in case of larger values of the yield stress. 

Solutions for the problem, using different computational techniques and different values of the yield 
stress, have been published by Begis m, (see also [29l Ch 6]), Sanchez [52], Mitsoulis and Zisis [53|, Vola et 
al. IH], Yu and Wachs [15|, Olshanskii m, De los Reyes and Gonzalez Andrade m, Zhang [5S], Glowinski 
and Wachs [S5], dos Santos et al. Syrakos et al. [S5], Aposporidis et al. [S5] and Muravieva [5D|. For 
specific features of the flow, such as vortex positions and intensities, we also refer to these works. 

From Figure]^ we observe that Algorithm FISTA* identifies the unyielded regions in agreement with the 
results published in the works cited above. The approximation computed with FISTA* including re-starts 
is overall similar. Nevertheless, the relatively large areas where the stress is very close to the yield stress 
make it difficult to detect where the stagnant flow region ends and where shearing begins. Overall, elements 
where \Th\ < Bi clearly dominate in these areas, which should, indeed, be classified as unyielded. 

Despite the identical stopping criterion in all cases. Algorithms ADMM and ISTA* clearly underestimate 
the regions occupied by unyielded fluid. While the approximation of the stress lies at least reasonably close 
to the yield stress in the blue areas, these two methods still fail to identify these as solid. 

Model Reduction with Adaptive Finite Elements. In past years, solutions on adaptive grids have already been 
successful at resolving the liquid-solid interface in fine detail, while reducing the substantial computational 
cost of simulations on uniform grids with the same fine resolution [^|6T1|^|^. Similarly, our objective is 
to achieve a resolution oi h = 1/128 in critical areas, while using a much coarser mesh with h = 1/16 where 
the residual is already comparatively small. 

For now, let us use the following ad hoc strategy: 

• Solve the optimisation problem with one of the four algorithms until convergence. 

(k) (k) 

• Determine the 60th percentile of the Frobenius norm of the residual IT’m/ ~ '7h \ triangles 

and refine those ^ 40% of all triangles with the largest residual. Further refinements of neighbouring 
triangles are required to avoid hanging nodes. 
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Figure 9: FrobeniuS norm of the stress |t| near the Bingham number Bi. Top left to bottom right: ADMM, ISTA*, FISTA* 
and FISTA* with re-starts. Values outside the range of the colourbar have been projected onto the upper and lower end points, 
respectively (Bi = 20, h = 1/128). 
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Figure 10: Adaptive finite elements for resolving the free boundary between yielded and unyielded regions. Stress magnitude 
|t^| computed by ALG2 (left) and FISTA* (right) for the same problem as in Figure]^ (Bi = 20). 


• Interpolate the converged solution linearly to the refined grid. 

In Figure we tackled the problem with Bi = 20 with ALG2 and FISTA* once again, this time on a grid 
that was only locally refined. Starting from the uniform mesh with h = 1/16, we cycled through the above 
refinement procedure three times. We conclude that the quality of both results is very much comparable to 
the one of the corresponding graphs in Figure]^ Nevertheless, it took about 65% (ALG2) or 61% (FISTA*) 
less computing time, respectively, until convergence was achieved. Additionally, the identification of the 
zero-flow region by ALG2 has even improved considerably. The upper stagnant zone still exhibits many 
coarse artefacts, though. 

Now that we can assume our refinement methodology to be validated for FISTA*, we apply it to the 
traditional challenge of predicting the yield surfaces when the yield stress is very large. We pick the two 
values Bi = 200 and Bi = 500 for which results have been published in the literature. We adaptively refine 


the initial homogeneous mesh with h = 1/16 five times. Our results are depicted in Figure 11 

Yu and Wachs [35] and Muravieva m have used ALG2 on a homogeneous grid with h = 1 /256 to solve 
these two problems. Their results deviate from each other as well as from ours, which are in close qualitative 
agreement with the publications of Mitsoulis and Zisis [S3] and Syrakos et al. [SS] . Since both of the latter 
works solve regularised approximations of the Bingham flow problem, fine geometric features like sharp tips 
that are visible in our results, have already been smoothed out in the problem formulations of these authors. 

Although this very basic approach to mesh adaptivity has proven to be effective, we anticipate even 
further improvements from more sophisticated, goal-oriented adaptive finite element methods like the DWR 
(dual weighted residual) method. We refer to the book of Suttmeier [^ for more details. 


/./. Conclusions 

We wish to emphasise that unlike second-order methods for instance, which require a Hessian at every 
iteration, the higher rate of convergence of the accelerated dual proximal gradient method FISTA* compared 
to the classical alternating direction method of multipliers ADMM / ALG2 comes at the very minimal cost 
of additionally evaluating a linear combination and storing one extra variable. Even though some open 
questions remain regarding the convergence of both FISTA* and ALG2, these seem to be of no practical 
relevance. For FISTA*, we have presented strategies related to the definition of the primal sequences, which 
allow us to generally prove an accelerated worst-case convergence rate of 0(l/fc), compared to only 0 {l/'/k) 
for classical algorithms. 
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Figure 11: Predicting the yield surface for flows at high yield stress values with FISTA* and adaptive finite elements. Top: 
|t^| and the corresponding grid when Bi = 200. Bottom: Bi = 500. 
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Furthermore, globally optimal values for free parameters that occur in the dual proximal gradient meth¬ 
ods can either be calculated a priori, or estimated numerically by backtracking. Moreover, the new dual 
FISTA method is very closely related to ALG2, in the sense that the subproblems that arise in both algo¬ 
rithms are either identical or even simpler for FISTA*. Any existing code based on an augmented Lagrangian 
method can therefore easily be modified to implement FISTA*. This is what leads us to our conclusion that 
the Algorithm FISTA* could be seen as a more efficient successor algorithm of ALG2 for solving genuinely 
nonsmooth formulations of viscoplastic flow problems accurately. 

We believe that at this stage it is still too early to express a recommendation towards either the simple 
FISTA* method or the variant with adaptive re-starting. Our small number of examples so far indicate 
that re-starting may not be be quite as effective in the context of Bingham flow as it is for other nonsmooth 
optimisation problems m- Further numerical studies are required to provide more guidance on this question. 

Our analysis already applies to a more general framework, including viscoplastic flow in three spatial 
dimensions. The extension to time-dependent flow problems follows naturally by first applying a suitable 
semi-discretisation scheme in time. If the inertial term is is discretised explicitly, then, at every time step, the 
methodology of this paper remains applicable except that an additional mass matrix arises in the solution 
of the Stokes problems. This strategy is completely analogous to solutions by other numerical methods, e.g. 

[nuBD]. 

We still see potential for improving the efficiency of FISTA* further, e.g. by employing preconditioning 
techniques (cf [65]) or inexact evaluations of the proximal map (cf [MlEl]). These concepts shall be our 
focus of further research on the topic. 
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